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A BLOCK ITERATIVE FINITE ELEMENT ALGORITHM FOR NUMERICAL SOLUTION 
OF THE STEADY-STATE, COMPRESSIBLE NAVIER-STOKES EQUATIONS 

By 

Charlie H. Cooke ^ 


SUMMARY 


An iterative method for niomerically solving the time independent 
Navier-Stokes equations for viscous compressible flows is presented. 
The method is based upon partial application of the Gauss-Seidel 
principle in block form to the systems of nonlinear algebraic equa- 
tions which arise in construction of finite element (Galerkin) models 
approximating solutions of fluid dynamic problems. The C®-cubic 
element on triangles is employed for function approximation. Compu- 
tational results for a free shear flow at Re = 1000 indicate signfi- 
cant achievement of economy in iterative convergence rate over finite 
element and finite difference models which employ the customary time 
dependent equations and asymptotic time marching procedure to steady 
solution. Numerical results are in excellent agreement with those 
obtained for the same test problem employing time marching finite 
element and finite difference solution techniques. 

INTRODUCTION 


Historically, most numerical methods for obtaining steady-state 
solutions of the nonlinear Navier-Stokes equations of fluid dynamics 
are time-dependent methods in which the steady solution is approached 
asymptotically by time marching procedures. Prevalence of such methods 
derives from the confidence that is placed in eventual convergence to 
the steady solution (for stable and consistent algorithms) , as well as 
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the ease and economy with which numerous finite difference methods may 
be formulated and implemented. However, in many situations the struc- 
ture of the problem or stability restrictions of the algorithm can 
lead to the requirement of many time steps to convergence. Hence, 
it is desirable to devise methods in which iterative convergence is 
achieved as rapidly as possible. 

In a series of investigations by Roache (refs. 1 and 2) several 
iterative methods have been devised for solving the incompressible 
2D steady Navier-Stokes equations in stream functions and vorticity. 

These finite difference techniques are neither time dependent nor even 
timelike in their iterations. Instead, the nature of the iterations 
is such that otherwise nonlinear equations approximating time independ- 
ent Navier-Stokes flows become linear in the independent variables 
over one iterative step. Recent advances in solving the resulting 
linear systems by direct methods (fast Poisson solvers and biharmonic 
solvers) are employed for economical equation solving. Improved 
iterative convergence rates have been experienced for the driven 
cavity problem and other low Reynolds number flows. 

On the other hand, applications of the finite element method in 
fluid dynamics problems governed by the unsteady Navier-Stokes equa- 
tions leaves one unconvinced such methods are at all competitive, 
in terms of economy of computer resources, with the standard finite 
difference techniques (refs. 3 and 4). A built-in cumbersomeness 
due to the general (gridwise) applicability of the method seems to 
indicate more complex program structure, more lengthy development and 
slower problem execution times are to be expected. Hence, in attempt- 
ing to make the finite element method more nearly competitive with 
the state-of-the-art finite difference algorithms which now exist 
in rather streamlined form, one might logically consider iterative 
solution procedures for the steady equations, with expectations that 

order of magnitude improvements in the number of iterations to con- 

• 

vergence could vastly improve the economy of the method. 

Early investigations of this type on the 2D steady transonic flow 
equations (ref. 5) are encouraging. A reported comparison of the Galerkin 
method (with quadratic rectangular elements) versus finite difference 
results for the flow around a circular-arc airfoil with imbedded shock 
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indicates finite element results which are superior to finite dif- 
ferences, while retaining a speed ratio of at least one order of 
magnitude. Here, of course, only one governing equation is involved, 
which greatly simplifies the numerical computation. 

Attempted solutions by finite elements of steady problems involv- 
ing the full Navier-Stokes equations have as yet been only sparsely 
reported. The case of viscous incompressible flow and laminar, steady, 
isothermal fluid motion in two dimensions has been investigated by 
Garling (ref. 3) , using quadratic function approximation with isopara- 
metric quadrilateral and triangular finite elements. The nonlinear 
algebraic equations resulting from the discretization process are 
solved with a Picard iteration of the form 

- 

where X is a vector of all density and momentum variables at the 
nodes of the discretization. The initial flow field employed the 
solution to the creeping flow problem (zero Reynolds number) and 
iterative convergence was achieved in only a few iterations, for a 
broad class of viscous flow problems and a significant range of 
Reynolds number (£ 10,000). The matrix A and vector f were 
assembled by the frontal solution technique, generalized to the case 
of nonsymmetric matrices (essentially Gaussian elimination without 
pivoting). [Herein lies the weakness of this formulation, since in 
the present investigation for the case of mixed subsonic-supersonic 
compressible flow at Re = 1000 we have experienced significant round 
off when Gaussian elimination with no pivoting (Grouts Method) was 
used. ] 

Laskaris (ref. 6) has developed a finite element numerical tech- 
nique whereby the steady state hydrodynamic equations for two-dimen- 
sional viscous compressible flows are solved, taking into full account 
the nonlinear convective terms, viscous terms, heat conduction terms, 
and variable fluid properties. The (Galerkin) method of weighted 
residuals is applied over distorted rectangular elements with cubic 
(Hermite, tensor product) function approximation. The resulting set 
of nonlinear algebraic equations for the nodal parameters are solved 
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by means of a multi-dimensional Newton-Raphson scheme. For a heat 
transfer problem in a diverging channel with plane walls, with around 
400 elements (roughly 2400 unknowns ) , convergence was achieved in 
5 to 8 iterations and approximately two hours CPU time on the GE-600 
computer. Direct Gaussian elimination and an out-of-core solver were 
used for matrix inversion at each iterative step. 

In the present investigation viscous compressible flows governed 
by the two-dimensional steady Navier-Stokes equations in primitive 
variables form are considered. The goal of the investigation is to 
determine whether the finite element method becomes a more feasible 
tool for fluids computations (for which either time asymptotic or 
steady governing equation formulation is applicable) when the steady- 
state governing equations are adopted. Solutions of the same physical 
problem by both time transient finite element and finite difference 
methods affords a ready basis for comparison of these diverse techniques. 

FLUID DYNAMICS MODEL OF A FREE SHEAR LAYER FLOW 

A computer code for numerical computation of 2D viscous fluid 
flows governed by the steady-state compressible Navier-Stokes equations 
in primitive variable form has been developed. Proof of concept for 
this finite element program is provided by the computation of the 
solution to a free shear flow generated by the parallel mixing of two 
supersonic jets, initially separated by a thin splitter plate. Solu- 
tion of the problem does not require the full Navier-Stokes equations, 
since fairly accurate results can be obtained using the quasi-parallel 
assumptions of parabolic boundary layer theory. However, the avail- 
ability of solutions generated by several computational methods (refs. 

4 and 7) affords a ready basis for evaluation of the finite element 
method. 


Flow Field Configuration • 

The flow field configuration of the test problem considered is 
shown in figure 1. The computational domain begins downstream from 
the base of the splitter plates. Numerical computations have been 
performed for flow at Reynolds number 1000. 
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Governing Equations 


Steady-> state flow is obtained through solution of the time inde- 
pendent Navier-Stokes equations. The assumption of constant total 
temperature (adiabatic mixing) and two-dimensional flow yields the 
following non-dimensional systems of governing equations (non-conserva- 
tive form) : 


Continuity 



+ u 


ap 

3x 


0 


( 1 ) 


y-momentum 



( 2 ) 


(3) 


Temperature relatio n 
T = 1 - u^ - v^ 


(4) 


Constitutive relationship 

Sutherland's viscosity law: y 


^3 2 


T + 198.6 1 
s 

T T + 198.^ 

- ® 


(5) 


5 



( 6 ) 


Perfect gas lawt P 


P 



T. 


In equation (5)/ u^T are dimensionless, although T and the 

9 

constant 198.6 (Sutherland's constant) are expressed in degrees 
Rankine. The variables used to non-dimensionalize eqs. (1) to (6) 
are presented in reference 7. 

Boundary Conditions 

Boundary conditions for the problem are shown schematically in 
figure 2. function specifications are given for all three variables 
on the inflow; symmetry conditions apply at the bottom, and on the 
top function specification is made for velocity, zero normal deriva- 
tive for density . On the outflow a computational boundary condition 
must be applied; this was chosen to be quadratic extrapolation. 

FINITE ELEMENT APPROXIMATION 


Our approximation to the fluid dynamics problem [equations (1) 
to (7) and boundary conditions of figure 2] is obtained by applying 
the classical Galerkin (or method of weighted residuals) in conjunc- 
tion with finite elements. The first step is to triangulate the 
computational domain ii with boundary T, and then consider piece- 
wise polynomial approximating (trial) functions on this grid. 


Trial Functions 

Function approximation in all independent variables is accom- 
plished by means of piecewise cubic trial functions on triangular 
elements. For a precise description of the element used, see refer- 
ence 8. For purposes of illustration, the trial functions for approxi- 
mating density variations are of the form 


P (x,y) 


N 

E Pj ♦j(x.y) 


(7) 
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Here N is the total number of nodes. The functions comprise 

a local and interpolating basis for functions of the form exhibited by 
eq. <7) which are continuous and piecewise cubic polynomial on n with 
sectionally continuous first partial derivatives , but which are infin- 
itely differentiable cubic on the interior of each triangle. 

The weights Pj are chosen by Galerkin's methods. Thus» the 
final approximating function satisfies all boundary and initial condi- 
tions at problem nodes / and approximately satisfies the governing 
equations over the domain. For each trial function (density and two 
velocity components) there are ten nodes per triangle; triple nodes 
at triangle vertices, and a single node at the centroid (see fig. 3). 
The parameters p. each associate with a distinct node^ and repre- 
sent approximations to function and first partial derivative values 

^p, "l^, at vertices, and function values alone at the centroid. 

The trial functions for velocity, defined similarly to those of 
density, are of the form 


N 

u(x,y) =* 53 Uj 4>j(x»y) 


N 

v(x,y) * I] v^ 4>-(x,y) 
J=1 ^ ^ 


( 8 ) 


Discretized Equations 

Consider a node at which a density independent variable Pj is 
not restricted by any boundary condition specification. The corres- 
ponding discretized finite element equation associated with Pj is 

determined by setting to zero the weighted residual obtained upon 

* 

multiplying eq. (1) by the basis function and integrating the 

result over fl. After shifting derivatives where possible onto the 
basis function (integrating by parts) and letting F»u,v be the 
vectors whose components consist of all nodal variables not restricted 
by any boundary condition, there results the system of equations: 
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(9) 


C(u,v)p » P(u,v) 

for determing the density of unknowns F* Here the typical element 
of C is specified by 


'JK 


II (“ 

Q 




J j^4»j4»j^(udx - udy)J (10) 


and 


^JL 


( 11 ) 


The sum defining Fj is over all nodes at which is known. The 

system of discretized equations for the velocity vectors u,v may be 
shown to have the functional form 


zz (u,v,p) 

1 zR(u,v,p) 


u 

Rz (u,V/p) 

1 

[ RR(u#v,‘p) 


V 


G (u,v,p) 
H (u,v,p) 


( 12 ) 


These equations are the result of multiplying eqs. (2) and (3) by 
and setting to zero the weighted residuftls obtained upon integratin«j 
the results over Q, employing integration by parts to shift (where 
second derivatives occur) higher derivatives onto the basis functions. 

The matrix elements of the u-momentuin matrices are specified by 
the equations: 



• (13) 
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( 14 ) 


^’‘JK 



- 2 + 


Ty" TnrJ“ 



»*K 

Sx 


dx 



‘j ^ 



» 


■// 

n 


P -5^ dA 


-/ 


P dy -£ (zzjj_ Uj_ + zBjj_ Vj_) . (15) 


The sum in eq. (15) is over all nodes at which u,v are bpeci- 
fied by boundary conditions. The v-momentum equations may be obtained 
from eqs. (13) to (15) by interchanging the variables x with y; 
u with V K.TiO reversing algebraic sign boundary integral terms. 
The indices OK in eqs. (10) co (15) range over all nodes at which 
P/U,v are unknown. 

To simplify accounting in the equation solving process, it is 
assumed that u„/V are either both known or both unknown at each 

Iv 1\ 

node. When one but not both is specified by boundary conditions# 
both are treated as unknown in the equation assembly process# and a 
corrected equation is inserted for the known component prior to equa- 
tion solution. This artifice produces momentum matrices character- 
ized by identical dimension# bandwidth# profile# and intraband distri- 
bution of zero and non-zero elements. Descriptors of any one matrix 
which must be stored then suffice for all. 

Numerical (Quadratures 

The appearance of nonlinearities in the integrands of eqs. (10) 
to (15) requires numerical quadratures. In order to maintain the 
degree of accuracy naturally achievable with cubic function approxima- 
tion# the analysis of Fix (ref. 9) which predicts a quadrature soheme 
exact for polynomials of total degree five is needed. These require- 
ments are conservatively met by the 16-point scheme (ref., 10) used by 
the finite element progreun# which possesses seventh order accuracy. 

A 7-point fifth order scheme was also tested and found to yield 
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approximately eight seconds per iterative step decrease in CPU time, 
with equivalent algorithm performance. 

A NONLINEAR BLOCK ITERATIVE GAUSS-SEIDEL SOLVER 

Consider now the problem of solving the systems (9) and (12) of 
nonlinear algebraic equations governing the approximating finite ele~ 
me/;t solution. Although more rapid convergence might be expected, it 
is clear that a full Newton iteration on these equations would be 
characterized by coupling of both continuity and momentum variables, 
implying excessive core requirements for in-core equation solving. 

As well, this method would lead to a very time consuming equation 
assembly process, since the Jacobian matrix would have rather compli- 
cated equations describing its elements. Ordinary nonlinear point 
SOR is not readily applicable* if a triangle by triangle assembly is 
desired; nor could convergence be readily guaranteed. 

However, a partial block Gauss-Seidel iteration* which is linear 
at each iterative step and which uncouples the density and momentxim 
solutions is readily designed. The equations proposed for this 
iteration are 


zz 

n 

RR 

n 


n+1 


u 


n+1 


n+1 


F 

n 


^n,n+l 

^n,n+l 


zR 

n 

Rz 

n 



(16) 

(17) 

(18) 


Here is used in assembling G,H, in order to update the momen- 

tum solution. One merit of this iteration is that only one system 
matrix at a time need be in memory during the equation solving. This 
greatly reduces core demands and avoids out-of-core solvers. 


* A full block Gauss-Seidel iteration cannot be applied, assuming one 
wishes for purposes of economy to simultaneously assemble both 
momentum eo'^ations. 
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In the present program eq. (16) is assembled first, with the 
matrix C in core. The resulting solution vector is input 

to the (simultaneous) assembly of eqs. (17) and (18) with RR in 
core and element matrices for zz on disk. Equation (18) is 
solved, then zz is assembled in core and eq. (17) is solved. 

The matrix inversions, for the supersonic free shear flow problem 
reported, were accomplished by Grout's method of LU decomposition. 


NUMERICAL RESULTS FOR A FREE SHEAR LAYER FLOW 


Case I. Time Independent Equations and Supersonic-Supersonic Inflow 

Steady- state results for computational solution of the supersonic 
free shear flow problem [eqs. (1) to (7) and boundary conditions of 
fig. 2] are now -lesented. The test case computecf employed Re = 1000. 
For a mesh confjit>;.ing of 225 elements (696 nodes per independent vari- 
able) iterative convergence was achieved with 15 iterations of eqs, 

(16) to (18) and 2079 seconds of CPU time on the CDC-6600 computer. 
Program core storage requirements were 162 Kg, For the 16-point 
quadrature scheme approximately 120 seconds per iterative step of 
CPU time were required, with approximately 8 seconds per step decrease 
when the 7-point scheme was employed. Total dollar costs for this 
computer run were $462, as determined by the computer systems' account- 
ing scheme. The convergence criterion used is 


Max 


“n 


< .001 


(19) 


where Af = f - f , and f is a function value of u, v, or 
n n+i n 

p, (Generally, the derivative values are less accurately modelled 
and lag function values in convergence.) 

For purposes of comparison these results may be considered in 
relation to those obtained for the same problem by various finite 
difference (ref. 7) and finite > lement approximations (ref. 4) applied 
to the time dependent Navier-atoke^ equations. Each method was initial- 
ized with the same starting flow field, and accuracy of the results 
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appears equivalent. Several indicators of computational efficiency 
for these methods are presented in table I. (The number of steps 
to convergence for the time transient finite element code is not 
the best obtainable since the maximum permissible step was not 
consistently applied, as for the finite difference runs. However, 
due to the expense of this method no attempts to rerun with maximum 
step were made. ) 

Figure 4 shows typical comparisons of steady- state density and 
velocity variations at stations Xj = 0.75 and X 2 = .175 in the flow 
(the strearawise extent of the computational domain is 0 5 x £ .225.), 
for steady finite elements versus time transient ADI finite differ- 
ences. Table II exhibits actual numerical differences between these 
computations. Table III presents percent differences with the ADI 
computations as base. Since the normal component of velocity is 
zero over much of the field, percent differences for this component 
are normalized with respect to the maximum value. 

Case II. Time Independent Equations and Subsonic-Supersonic Inflow 

The steady finite element code has also been applied to a free 
shear flow problem resulting from the mixing of a subsonic and a 
supersonic flow. For a description of this problem, refer to figure 
1 with Mj = 0.11, M 2 = 3.00, and the boundary conditions of figure 5. 

Catastrophic failure of the method for this case resulted from 
the equation solving, arising from very ill-conditioned matrices and 
Gaussian elimination without pivoting. For example, a standard de- 
bugging procedure is to apply the code to a constant flow problem 
(say, p = .07625, u = .8018, v = 0; or v = .8018, u = 0, for y- 
direction flow) to see if this flow reproduces itself. With large 
meshes (around 300 triangles) thought necessary for the subsonic- 
supersonic case the flow did not reproduce in the u-component, for 
y-direction flow. It was determined that the system matrices pos- 
sessed determinants of order of magnitude 


det(zz) = lO"^®®^ 
det(RR) = 10“®®3! 
det(C) = 10-1775 


( 20 ) 
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This seems to indicate ill-conditioning of all, but with a much more 
severe occurrence for the u-momentum equations (observe that the ratio 
of det(zz) to det(RR) is 0(10"^®). 

As a further test for the ill-conditioning, the equation 

zzu = ¥ (21) 

was, by proper choice of T, set to have a solution all of whose 
components assumed the value of 1. This system was then solved by 
Grout's method (standard finite element solver) and by a band solver 
which used partial pivoting. Partial pivoting gave the best results, 
but each solution was characterized by numerous values with no signifi- 
cant figures of accuracy. Neither solver produced diagnostics peculiar 
to a singular matrix, although the eqs. (20) indicate near singularity. 
Best results, although still unacceptable, were obtained from Grout's 
method with double precision inner product accumulation on the decompo- 
sition, with matrix rows scaled so as to have unity diagonal elements. 

As a remark on how nearly singular a matrix can be and the systems 
(16) to (18) still be solvable by ordinary techniques, the computa- 
tional success for the supersonic-supersonic case prevailed in spite 
of the characteristics of near singularity indicated by 

det(zz) ^det(RR) q/10"1700 (22) 

To make the steady finite element formulation a generally appli- 
cable tool, these ill-conditioning problems must be overcome, and the 
feasibility of the resulting method evaluated, 

Gase III, Time Dependent Equations and Subsonic-Supersonic Inflow 

During the present work period the finite element code developed 
previously (ref. 4) for numerical solution of the time dependent* 
Navier-Stokes equations has also been applied to the free shear flow 
resulting from the mixing of subsonic and supersonic streams. 
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Figures 6 shows comparison of the finite element results (after 
600 steps# At ^ .01) and steady state results for the central dif- 
ference ADI code. Data output from the two codes compare well with 
the exception of the normal component of velocity# which exhibits a 
maximum difference of 4 to 6 percent (relative percent difference) 
near the top right corner of the flow field. Here the finite element 
computation had not completely converged. This local slow convergence 
is attributed to boundary condition errors on the top; the finite 
element flow domain was obtained by truncating a region 20 finite 
difference mesh increments in width from the top of the ADI domain. 
Boundary conditions were then supplied by ADI steady results at the 
finite element domain. Initially# a significant bulge in finite 
element data occurred in this region# with the rest of the field con- 
verged to steady state. After the ADI code had been run further in 
time with a more solid convergence check instituted# the top boundary 
conditions then supplied resulted in a rapid and significant decrease 
in the bulge occurring in finite element output# to the level now 
indicated. It is conjectured mesh refinement of the ADI domain would 
be required in order to produce boundary condition data sufficiently 
accurate to produce totally satisfactory global convergence of the 
finite element calculation. Further support for this conclusion is 
furnished by figure 8 of reference 7# which shows disagreement in the 
inviscid region between various finite difference solutions# implying 
inaccuracy near the. top of the finite element domain. 

The time transient finite element code was also applied to a 
high Reynolds number (Re = 80# 625) supersonic-supersonic mixing problem 
(see figs. 1 and 2). At a time step of .01 catastrophic failure 
emerged rapidly (negative temperatures in 30 steps) . At a time step 
of .001 the computation was running smoothly at 150 steps# with the 
usual convergence criteria satisfied. However# these convergence 
criteria appeared satisfied through most of the run; it was conc]^uded 

the initial flow was also steady state or else the time span insuffi- 
cient for significant changes to develop. From this result it appears 
high Reynolds number flows without shocks could probably be calcu- 
lated with this code. 
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CONCLUSIONS 


A finite element code for numerical solution of fluid flow prob- 
lems characterized by the two-dimensional time independent Navier- 
Stokes equations has been developed. Proof of concept was provided 
by the calculation of the primitive flow variables for a free shear 
flow problem. Excellent numerical results were obtained in compari- 
, so.'* to ADI and various other finite difference methods. 

For the supersonic-supersonic free shear layer problem, order of 
magnitude improvement in iterative convergence rate (15 steps compared 
to over 100 steps) was achieved, in comparison to the previously 
developed time transient finite element code (ref. 4), with some 
improvement in storage (236 Kg down to 162 Kg) over the most feasible 
version of this code. Moreover, reduction in number of equation terms 
due to the steady form of the governing equation, as well as the 
diverse natures of the two numerical processes, produced reductions 
in CPU time (154 seconds/iterative step down to 120 seconds/iterative 
step) and 0/S calls (7166/step down to 1025/step). From all these 
factors there resulted a reduction in machine total dollar cost (as 
calculated by the CDC-6600 operating system's accounting routine) 
for obtaining the converged solution on the order of 20 to 1, not to 
mention the significant reduction in man-hours necessary for process- 
ing the multiplicity of computer runs required by the time transient 
code. Thus we may readily conclude that the steady-state finite 
element approach is far more efficient than the tii.ie transient 

% 

formulation. 

On the other hand, at this stage of the investigation, there 
apparently exists the drawback of a less general applicability of 
the steady formulation; the weakness exhibited (ill-conditioning) 
for the 300-element mesh and the subsonic-supersonic flow problem. 

For broad general applicability of the code, the problem of ill- ^ 
conditioning has yet to be overcome, and the resulting feasibility 
of the method then evaluated. This problem appears to be related to 
the mesh size; certainly the system matrices become more nearly 
singular as the number of nodes increases, and roundoff effects 
have further room to propagate. 
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Moreover, even with the success of the steady finite element code 
in evidence when compared with the time transient finite element code, 
at this state of development this method is still not quite competi- 
tive with the better finite difference techniques. However, the gap 
has been significantly narrowed to the point where the finite element 
method can almost be considered a feasible alternative. 

As regards the time transient finite element code, it appears to 
have provided a reliable computational tool, for all test problems to 
which it has been applied, at the expense of much too heavy demands 
on computer resources. 


VI 


EPILOGUE 

Several factors contribute unnecessarily to a broadening of the 
gap between finite element and finite difference results. For exampler 
it has been determined that equation assembly time in comparison to 
equation solving time per step is highly unbalanced on the side of 
equation assembly time. This imbalance could be lessened several ways: 

1. Triangular elements were employed solely for general appli- 
cability of the finished code to other than rectangular regions. How- 
ever, the code is being measured for competitiveness using a test 
problem whose domain is a rectangular region. For such a case, the 
(tensor product) Hermite cubic shape functions on rectangular elements 
would lead to more sufficient equation assembly time in at least two 
ways — fewer* (one-half) as many elements to process for the same mesh 
accuracy; and the existence of more efficient quadrature schemes for 
rectangles than for triangles (the seventh order accuracy afforded by 
the 16-point triangle scheme is afforded for a rectangle by an 8- or 
9-point scheme) . Finally, significant finite element results thus far 
(refs. 3, 5, and 6) have been with rectangular elements. 

2. A scheme originally intended to improve the efficiency of the 
method actually degrades it, on the momentum equation assembly — for 
example, integrals like 

puv dA (23) 

can be evaluated each iterative step from p,u,v values at the quadrature 
points (which must be computed anyway to determine viscosity) by apply- 
ing the quadrature scheme directly to eq. (23), or the integrals 



* Two triangles composing one rectangle. 
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(24) 


may be computed once and stored, for future computations at each 
step of the form 


// puv dA -L ")) 


(25) 


For an integrand with an independent variable product of degree 
higher than two the second scheme is less efficient than the first, 
assviming the 16-point quadrature scheme. Consequently, in terms of 
number of multiplications necessary the momentum equations assembly 
suffers from inefficiency, possibly a significant amount. 

3. It is conjectured that a linear element code would be more 
efficient on the equation assembly. Here only a 1-point quadrature 
scheme is necessary for the accuracy needed (ref. 9), a great simpli- 
fication. However, more (triangles) elements would be required. One 
could only guess whether the iteration scheme would converge as well 
and what effect alternate boundary condition implementations would 
cause. 
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TABLE 1- CONVERGENCE REQUIREMENTS FOR VARIOUS COMPUTATIONAL METHODS 
(Re - 1000, SUPERSONIC-SUPERSONIC FREE SHEAR LAYER FLOW) 
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^ Not best number of steps , since maximum allowable step was not consistently used 























TABLE K- NUMERICAL DIFFERENCE BETWEEN FINITE ELEMENT AND ADI STEADY STATE RESULTS 

ISUPERSOmC-SUPERSONIC FLOW) 
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TABLE m.- PERCENT DIFFERENCE BETWEEN FINITE ELEMENT AND ADI STEADY STATE RESULTS 

(SUPERSONIC-SUPERSONIC FLOW) 
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- Flow field configuration, supersonic-supersonic shear layer problem. 
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Figure 2.- Boundary conditions for supersonic-supersonic flow. 
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Figure 3. - Node numbering scheme for the C® cubic element. 
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Figure 5. - Boundary conditions for mixed subsonic-supersonic flow. 


o Finite Elements 



(a) X • 0. 15 

♦ * 
Figure d.- Comparison of time transient Finite Element and ADI steady state results; 
mixed subsonic-supersonic flow. 



o Finite Elements 




>. 



(b)x- a35 

Figure 6. - Concluded. 
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